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We investigate the quench of half-filled ID and 2D fermionic Hubbard models to models without 
Coulomb interaction. Since the time propagation is gaussian we can use a variety of time-dependent 
quantum Monte Carlo methods to tackle this problem without generating a dynamical sign problem. 
Using a continuous time quantum Monte Carlo method (CTQMC) we achieve a system size of 
128 sites in ID, and using a Blankenbecler-Scalapino-Sugar (BSS) type algorithm we were able 
to simulate 20 x 20 square lattices. Applying these methods to study the dynamics after the 
quench, we observe that the final state of the system can be reasonably well described by a thermal 
single-particle density matrix that takes the initial single particle conservation laws into account. 
The characteristic decay towards this limit is found to be oscillatory with an additional power 
law decay that depends on the dimensionality. This numerically exact result is shown to compare 
favorable to mean-field approximations as well as to perturbation theory. Furthermore we observe 
the information propagation in the ID-case in the charge charge and spin spin correlations and find 
that it is linear with a velocity of roughly v ~ 4 in units of the hopping amplitude. 
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I. INTRODUCTION 

The realization of various solid state Hamiltonians us- 
ing ultra-cold atom experiments has sparked great inter- 
est in their out-of-equilibrium physics. Ongoing experi- 
mental work has achieved great control in the preparation 
of those systems and experimental physicists have started 
studying the dynamics of these systems when they e.g. 
change parameters of their trapping devices [TJ [2] . In ad- 
dition, pump and probe femtosecond spectroscopy per- 
mits the study of electron relaxation dynamics [3l 0] . In 
this article we propose quantum Monte Carlo methods 
(QMC) that allow the numerical real-time evolution of 
quantum systems which are prepared in an initial ther- 
mal state that contains arbitrary correlations. The initial 
thermal state is computed within the QMC algorithms, 
which arc formulated in terms of a path integral on the 
Matsubara-Keldysh contour. Of special importance is 
the fact that these methods allow for quenches to ar- 
bitrary non-interacting models without introducing an 
additional dynamical sign problem. In this paper we 
apply these methods to 1- and 2-dimensional half-filled 
Hubbard models prepared in a thermal initial state and 
quench them to a Hamiltonian without the Hubbard in- 
teraction, that is U(t > 0) = 0. Sotiriadis ;5. called this 
situation with a thermal initial density matrix a thermal 
quantum quench in contrast to the pure quantum quench 
of a pure initial state. 

In this setup, we can ask a number of questions con- 
cerning the evolution of the system after the quench. 

• Does the system evolve to a new steady-state? 

• How does an isolated system approach a possible 
new equilibrium? 

• What is the nature of this state? 



• Does the system retain memory of the initial state? 

In ID we use an extension of Rubtsov's CTQMC method 
[HE], and in 2D an extended BSS type algorithm [8]. A 
complementary model to ours, where, starting from the 
free electron limit, the Hubbard interaction was switched 
on at t — 0, was studied theoretically in [5] and nu- 
merically by Kollar et. al. [10] using dynamical mean- 
field theory (DMFT). Kollar et. al. found a critical 
value of Uc ~ 3.3 where the characteristic oscillations 
in the double occupancy and the fermi surface discon- 
tinuity seem to have been suppressed. They call this a 
dynamical phase transition as they observe a very fast 
thermalization in this regime. More general results for 
the quench dynamics of a quantum system in arbitrary 
dimensions have been presented by Moeckel and Kchrcin 
[llj . Manmana et. al. [H] studied a similar problem 
as ours using time-dependent density matrix renormal- 
ization group (DMRG) techniques but in contrast to our 
spinful electrons they considered the case of spin-less elec- 
trons. They quenched from interaction parameters lying 
in the metallic or insulating regime to specific values of 
final Hubbard U's, where they also crossed those phases. 
They found that the information propagation in their 
system, as observed in their density density correlation 
functions, happens only with a finite velocity that de- 
pends on the final interaction and not instantaneously. 
We also observe this finite velocity of propagation in the 
charge charge correlation functions for spinful fermions, 
which gives rise to the notion of a light cone like evo- 
lution of the information propagation. For a number of 
models this finite velocity of the propagation of infor- 
mation is known as the Lieb-Robinson bound and was 
first discovered for quantum spin systems [13] . Lieb and 
Robinson proved that in their system only exponentially 
small corrections exist outside this light cone. Over the 



years these theorems got extended to more systems up to 
arbitrary harmonic systems on general lattices with lo- 
cal dynamics (see Ref. [H] and references therein) . This 
light cone like structure is a direct consequence of the 
locality of the dynamics. 

The structure of the article and our main results are 
the following. We first define the Hamiltonian and its 
symmetries in |Sec. II| and then we carry out mean-field 
and perturbative calculations to gain insight into the 
physics at hand. This is summarized in |Sec. IlT| |Sec. iV| 
describes the two QMC algorithms used (CTQMC and 
BSS) to study the physics of the quench, in a self- 
contained manner. Our exact numerical results are pre- 
sented in |Sec. V| and compare favorably with the analytic 
calculations. We find that in ID and 2D local quanti- 
ties equilibrate to values that can be reasonably well de- 
scribed by an effective single particle density matrix that 
respects the particle densities n& of the initial thermal 
density matrix. The approach to the equilibrium follows 
a dimension dependent power law. For single particle 
quantities such as Green functions this power law follows 
that of a diffusion process, t~ D / 2 , where t is the time. 
Finally we take a look at the information propagation of 
correlations in real space through the system. We show 
that a light cone like structure exists, beyond which the 
propagation of the information of the correlation is ex- 
ponentially suppressed. Details of the calculations are 
presented in appendices. 

II. THE MODEL AND ITS SYMMETRIES 

Using the well-known fermionic operators where cj CT 
creates an electron and Ci a annihilates an electron we 
can define the Hubbard model: 

H = ~ H %4 c jv + u J2( n n - - \) • (!) 
- . ' * ' 

= :H a =:Hu 

Here = c ia Ci a . For our simulations, we restrict our- 
selves to the case of nearest-neighbor hopping on hyper- 
cubic lattices with the lattice constant set to unity. 

{f if i, 7 are nearest neighbors 
otherwise 

and we restrict ourselves to half-band filling. With this 
choice of parameters, the negative sign problem does not 
plague the evaluation of the initial density matrix. Defin- 
ing the hopping matrix in that way we set the energy unit 
to the amplitude of the hopping. This frees the letter t 
and enables us to use it for denoting the real time t. The 
above defines our unnormalized initial density matrix: 

p{t = 0) = e-< m (3) 

describing a Mott insulating state at inverse temperature 
j3. At t — we switch off the Hubbard interaction, that 
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FIG. 1. The momentum resolved rj - pairing correlation func- 
tion as a function of t for an initial Hubbard interaction of 
U = 2 on a L — 128 site chain. One can clearly see the hori- 
zontal black line in the middle of the figure which is present 
for all times. This corresponds to the fact that this special 
value of q — n is a conserved quantity. 

is Hjj(t > 0) = 0, such that the unitary time evolution 
of the system is given by 

U{t,t')=e-' lH ^ t ~ t '\ (4) 

As the evolution of the system is carried out with the hop- 
ping Hamiltonian Hq, there are a number of conserved 
quantities. In particular the fc-space resolved particle 
density and all related quantities, such as the kinetic 
energy, are conserved since [rika,Ho] = 0. 

In addition, at the particle-hole symmetric point, 
where 

e{k) = -e{k-Q) (5) 

holds, 77-pairing 

k 

is a conserved quantity, since 

r)l(t) = Y e lt « fe > £ (-*+Q)) c t c f - - (7) 

k 

D 

with e(k) = — 2 cos(fci) in D dimensions and 

i=i 

D 

Q = 7T e*j where e*j denotes cartesian unit vectors. 

i=i 

Therefore 77-pairing provides a non-trivial test of the 
QMC methods. |Fig. 1| shows the result of a Monte Carlo 
run. The black line that is present for all times, in the 
middle of the figure corresponds to the particle-hole sym- 
metric point in the 77-pairing. Its existence provides a 
non-trivial test for the QMC methods that we outline in 
ISec. IVI 
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III. MEAN-FIELD ANALYSIS AND 
PERTURBATION THEORY 

A. Mean-field analysis 

We carried out a mean-field approximation based on an 
anti-ferromagnetic decomposition of the Hubbard Hamil- 
tonian. In ID, the unit- vector ai = (1), and in 2D, the 

basis ai = (1,0) T and a 2 = (0, 1) T span the lattice. 
To have the possibility of an anti-ferromagnetic order- 
ing in the mean-field approximation we define the anti- 
ferromagnetic unit cell in ID simply as twice as large 
as A\ — (2) and in 2D we define the unit-vectors as 
Ai = (1, 1) T and A 2 = (1,-1) T . The unit cell now con- 
tains two orbitals and we label one of them with Cj and 
the other with di. We define the Fourier transform of 
these operators as 



cr 



e ikR c k 



e lkR d k 



(8) 



where R denotes lattice vectors labeling anti- 
ferromagnetic unit cells. With the definition of the 
anti-ferromagnetic unit cell and the labeling of its elec- 
trons we can now proceed to describe the calculation in 
a unified way. Additionally we introduce n d Ra = d Ra dncr 
as well as n Ra = c Ra CR a as shorthand notation. We 
rewrite the Hubbard interaction as 



H, 



-U 



L R\ 



(9) 



In this mean-field approximation we introduce the mean- 
field order parameter m z : 



m z = {n c m - n c Ri ) 
-m z = (n d m - n d Rl ) 

such that the mean-field Hubbard interaction reads 
Um z \ ~\ 



ttMF 



-y- 1^ [ n m - n Ri) 

R 



[n R ^ n d R ^) 



(10) 



(11) 



Carrying out the Fourier transform on the kinetic energy 
part and introducing the quantity ~j ka = (c ka ,d kcr ) T we 
can rewrite the total Hamiltonian in a matrix form like 



ka 



-Act Z k \ ( c ka 
Z k Act J \ d ko 



(12) 



= :H (k,a) 



Here Z{k) contains the information about the lat- 
tice structure, Z(k) = — (l + e~ lfcj4 ) in ID and 
Z(k) = -(l + e- ikAl ) (l + e- tkA *) in 2D. A = rnJL 
and the band structure is determined by E(k) — 



ztyj A 2 + \Z(k)\ 2 . This, so far, is well known from ther- 
modynamics. Next we need to introduce real-time dy- 
namics into this setting. The Hamiltonian responsible for 
the time evolution contains no explicit time-dependence, 
therefore the evolution of 7 is given by: 



e iHot lka e 



iH Q t 



(13) 



To obtain a Heisenberg equation of motion, we derive 
this equation with respect to t. Taking into account the 
structure of 7 as well as the fermionic commutation rules, 
we get 



dt 



ik.At) = -i(H a (k,*hl <T (t)) 



This is solved by 



7fc«r(t) = e- iff °^ ff ) 4 7fc CT - 



(14) 



(15) 



The time dependence of the magnetization is then given 
by 



ka 



(16) 



This quantity is plotted in Fig. 2 (a) as a function of di- 
mension D and at T = 0. From the log-log plot (Fig. 2 



(d)) we clearly see that the maxima of the oscillations 
can be fitted by functions that decay as a power-law. To 
get rid of certain numerical artifacts in the log -log-plo t of 
the absolute value we also plot m z (t) ■ t D l 2 in Fig. 2 (b) . 
We see that this quantity quickly approaches a simple 
sine wave oscillation pattern with constant amplitude. 
In Fig. 2 (c) we show the Fourier transform of this quan- 

4 is present 
. De- 



tity. In ID only a single frequency of uj 
whereas 2D oscillates with a frequency of ui 
pending on the dimensionality we see an odd-even effect. 
Even dimensions have the same base frequency as in 2D 
and odd dimensions have the same base frequency as ID 
has. Additionally, higher harmonics show up in D = 3 
and D = 4 with a spacing of Aw « 8. Judging from this 
plot we propose that the long-time behaviour of the en- 
velope of this decay is connected with the dimensionality 
of the system like 



\m z (t)\ oc t 



(17) 



In the thermodynamic limit we can perform a more de- 
tailed analysis of the behaviour of the ID magnetization. 
After some calculation we get 



m z (t) 



2A 



dx- 



cos(2£.t) 



V4 - x 2 y/x 2 + A 2 



(18) 



The limit A — ¥ 00 is exactly solvable and is a representa- 
tion of the Bessel function Jo, therefore m z (t) — Jo(4t). 
The leading order asymptotic behaviour of the ID mag- 
netization with respect to t is given by 



m z (A,t) = (l + —)-iJ (4t) 



(19) 
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Analysis of mean-field magnetization 
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FIG. 2. The time-dependent decay of the magnetization of 
the mean-fiel d solution with T — and U = 2. (a) shows the 
behaviour of Eq. (16) in dimensions from 1 to 4. The equa- 



tions for the magnetization were solved self-consistently and 
afterwards the time-propagation was calculated. The colors 
are consistent through all four plots. Now in (b) every data 

set is multiplied by the expected t~? behaviour. We see that 
all data sets approach sines with constant amplitude, thereby 
confirming our conjecture |Eq. (17)] without introducing the 
numerical artifacts of the log-log-plot (d). (c) is the Fourier 
transform of the data sets in (b). Finally (d) is the double 
logarithmic plot of the absolute value of the data sets in (a). 
The kinks in the data for ID and 3D are due to the fact that 
we have data with a very fine time-resolution. The oscillations 
in these graphs would extend all the way to zero in the plot, 
thereby effectively coloring the plot in black. Nevertheless the 
power-law decay of the envelopes is clearly visible. 



therefore confirming our hypothesis in ID. In Sec. A 2 
we generalized this to arbitrary dimensions and give an 
asymptotic expansion of m z (t,A) with respect to A. A 
key observation is that in dimension D the limit A — > oo 
is exactly solvable 

m D (t) := m z , D (*, A -> oo) oc Jj?(4*), (21) 

therefore the asymptotic behaviour to lowest order in 
time of the magnetization is 



m D (t) (xfT sm. D (it + -) 



(22) 



which fits well to our numerical observations at a finite 
A. The observed frequencies can be explained by us- 
ing power reduction formulas for trigonometric functions 
|15j . In the following 0Ju,k — 4(D — 2k) denotes the fc'th 
frequency. For even dimension we have 



sin D (4t) = 



1 



'—k 



k=0 



cos(wD :fe i) 



(23) 

Therefore in even dimensions the smallest observed fre- 
quency is 2wo where a>o is some base frequency(in our 
case we have u>q = 4) and the largest observed frequency 
is Dloq. To this observation fits a simpler result derived 
by using a fiat band of bandwidth 2w in |Sec. A3] The 
behaviour in the limit of A — > oo is 



m c z onBt (w,t) oc 



sin(2it)i) 
2wt 



(24) 



Hence the decay of the system with a constant density of 
states is similar to a 2D Hubbard system with bandwidth 
2w = 8. This result is consistent with the t~?~ law, since 
a constant density of states is realized by free electrons 
in a 2D continuum. Here it is obvious that the frequency 
of the oscillations depends on the bandwidth. In odd 
dimensions the powers of the sine are given by 



sm 



D 



(4i) 



1 



2 D 



-l A^i 

k=0 



sin(w Di fet). (25) 



In odd dimension we conclude that the lowest observable 
frequency is lvq and that the largest frequency is again 
Du Q . We see that in the large-i regime the frequencies 
are given by u>k = (D — 2k)ujQ. Therefore it is clear 
that the difference between two frequencies is Aw = 2cjq 
which gives in our case the observed Auj — 8. By noting 
that the squared magnetization is related to the double 
occupancy by 



N 



(26) 



The calculation is outlined in lScc. ATI The well-known 
asymptotic behaviour for large t of Jq is 



J (4t) = sin(4£ 



4' 



1 1 

/2tt y/i 



+ o(t- 3 / 2 ) 



(20) 



we expect to see similar behaviour in accessible correla- 
tion functions in our QMC data. We conclude this section 
by noting that for a fixed linear dimension size effects set 
in at the same time. This is consistent with the notion 
of a dimension independent velocity for the propagation 
of information. 
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B. Perturbation theory 



Comparing Monte Carlo data and perturbation theory 



Here we consider a first order expansion in the strength 
of the interaction in the thermal Hamiltonian. This is 
sufficient since no interaction is present in the real-time 
evolution. Additionally, the expansion has the advan- 
tage that we can analytically perform the thermodynamic 
limit. We consider the spin spin correlation function 

S(Ri,t) = oo'{n ia (t)n 0cr ,(t)) 

re ifli9 <4 CT (t)c fc _ g(7 (t) C t (t)^+^(t)) 



era' 



(Tcr'kpq 



S(q,t). 



aa' ' kpq 



(27) 



The usual perturbative expansion of p gives 



(S aa >(k,p, q,t)) = ((1 - / dTHu(T))S arT '(k,p,q,t)} . 



(28) 

This result is the same as obtained by an expansion 
of the full Keldysh evolution operator Sc- It is obvi- 
ous that all contributions stem from correlation func- 
tions that mix real-time and imaginary-time. Therefore 
a solution of this problem using a plain Keldysh method 
along the real-time contour is not possible. After Wick- 
decomposing this expression and collecting the remain- 
ing terms, we get for the spin spin correlation function 
S(q,t): 



'2 



2U 



JdrJ2 G< +9 (r, t)G>(t, r)G<_ q (T, t)G> (t, r). 

(29) 



kp 



To interpret the dimensional dependence in QMC sim- 
ulations we will take a closer look at S(w,t), since it 
is related to the magnetization. Neglecting the time- 
independent zeroth order contribution we get: 



S(n,t) = const. - 2U J dr£(t, r)£(i, r 
o 



deg(e)(f(pe) - 1)V<*+-) 



(30) 



(31) 



In the last line of Eq. (31) we performed the thermo- 




FIG. 3. The spin spin correlation function S(q = (iv, tt) t , t) of 
a 2D Monte Carlo simulation of a 20 x 20 lattice at U = 8 and 
/3 = 2. 5. We exp ect this to be in the high-temperature regime 
where Eq. (33)] is valid. The dashed black line is oc Jg(4f). 
The offset is taken from the large time behaviour of the QMC 
data and the amplitude was taken from the value at t = 0. 
The inset shows a magnified view of the region below the inset 
from t = 1 to t = 4.5. We see that the approximation works 
almost flawlessly in that regime. The deviation for t — > 5 has 
its root in boundary effects that set in for approaching that 
point in time. 



/(e) denotes the usual Fermi function. Performing the 
r-integral in S(t) and rearranging terms we get: 



oo oo 

S^tt, t) cx J de k J de 



g(e P )ff(e fc ) 2zt(£p+efc) , 
£fc + £p 



x [f(l3, ~e k )f 2 (l3, -e p ) - / 2 (/?, e,)/ 2 (/3, e p )] . 

(32) 

The large bracket that contains all Fermi functions has 
an expansion in j3 as § • (e p + e k ) + Q(/3 3 ). Therefore, to 



first order the denominator in Eq. (32) is canceled and 
the two remaining integrals decouple. Then we have 



S(ir, t) cx 



deg(e)e 



2id 



(33) 



Specializing to the density of states of a ID chain 
g lD (e) = ^Eff we get S(ir,t) cx J 2 (4i) with J the 
Bessel function of the hrst kind. Since the spin spin corre- 
lation in the high temperature limit is just the mean-held 
magnetization squared we can deduce from the general 
result in Sec. A 2 that the lea ding ord er decay of S(ir, t) is 



like S(ir,t) cx t~° sin 2iJ (4t). Fig. 3 shows a comparison 



dynamic limit and introduced the density of states g(e) 



of this theoretically predicted behaviour in this approxi- 
mation to an exact Monte Carlo run in 2D. We note that 
the large /? limit gives the same leading order behaviour. 
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IV. DESCRIPTION OF THE QMC 
ALGORITHMS 



The weak-coupling CTQMC approach for 
real-time evolution 



The physics of the Hubbard model is usually not rea- 
sonably well described by the simple approximations of 
the previous section. Especially the mean-field descrip- 
tion of the ID Hubbard model is usually just plain wrong 
since it fails to describe its low energy Luttinger liquid 
physics. Therefore we have to verify and extend our ap- 
proximate results using unbiased numerical QMC meth- 
ods. Starting from the average of an operator On{t) in 
the Heisenberg picture, we have the time-dependent av- 
erage 



0(t) = Tr (pO H (t)) 

= Tr(pU(0,t)O s (t)U(t,0)) 



(34) 



with the density matrix p of the Hubbard model 
Eq. (1) the time evolution operator U(t, t') and the 



possibly explicitly time-dependent operator Os(t) in the 
Schrodinger picture. This equation forms the basis of 
the stochastic methods outlined in this section. Insert- 
ing the identity 1 gives the partition function which 
allows the construction of the Markov chain Monte 
Carlo method in |Sec. IV A") Restricting the real-time 
dependence in U{t,t') to single particle Hamiltonians 
enables the treatment of time-dependent problems by 
using BSS type algorithms as outlined in |Sec. IV C| 



Starting from the partition function we derive a weak- 
coupling CTQMC method on the full contour similar as 
done in Refs. [H [16] on the imaginary contour. The 
partition function Z on the Keldysh contour is given by 



Z 
Z~ 



=£*= 

71 = 



dz r ■ ■ dz n (TcH u (z 1 )H u (z 2 ).-H u (z n )) Q 



c 



(35) 

where Tp orders the c ontour times along the Keldysh 
contour C (depicted in Fig. 4), Zq — Tre~^ H ° and the 
Hubbard interaction Hjj(z) in the interaction picture. 
Similar as in Ref. [7] we introduce an additional Ising 
spin s 1 into Hjj, hence 



Hr 



= t££( 



"m -1,-^)^-1+^6) (36) 



as well as the new parameter 5. From thermodynamic 
QMC it is known that 6 can be used to reduce the sign- 
problem of the simulation and we can confirm that it can 
still be used to tune the sign in the fully time-dependent 
setting. The usual choice for ID Hubbard models to 
eliminate the sign-problem is S = | + + . Introducing 
that into the general expansion for the partition function 
Eq. (35) gives 



Z 



2 / 



£ 

ra=0 



dzi £ " ' / dZn £ IJ(Tc(«ii,o-(2i) - cVs 1 ) ■ ■ • (n initr (zn) - a<r, S "))o 

ius 1 c 



(37) 



l n ,S n (7 



where we introduced a a „• = \ +<rs l 5 and made use of the where 



fact that for S^-conserving problems the weight splits up 
in an f-part and a J,-part. We can compactify Eq. (37) by 
introducing configurations. A configuration C' n consists 
of Hubbard vertices Vj — [ij, Zj, s- 7 ] with their Ising spin 
s- 7 , hence 



det(M ff (C„)) = 
G° til {z\,zi) - a a 

^i„,ii Zl ) 



^ii,i„ ( z li z n) 
ij, ( z nj Zn) &o 



The entries of M a (C„) are given by the free Green's func- 
With that concept we can introduce the sum over the Hon 



C„ = {[ii,zi,s 1 ],...,[i n ,z„,s™]} 

concept 
configuration space 

£ = £^ / dz i£--- / " dz " £ 

C„ n=0 ' i, i 1>s i q »„.a" 



(38) 



(41) 



(39) 



(T c cl(zi)c k (zk))o ~ S lk a a ,. 



(42) 



For the Monte Carlo evaluation of the contour integrals 
Eq. (39)| we have to transform them to linear integrals. 



Using this notati on and a pplying Wick's theorem to the To ac hi ev e that we need to specify the parametrization 
thermal average, |Ec~ (37)| can be rewritten as a sum over f tue contour. An obvious linear one is 
determinants 



Z 



£ 



V 



l[det(M a (C n )) (40) 



z(s) 



b exp 



s G [0, t exp ] 



(43) 



l(^S ^texp) S £ exp •> exp 



7 



C 



to t+ 
T 



to-iP 



FIG. 4. The full contour C that enables us to cover imaginary- 
time evolution and real-time evolution on a common footing. 
t~ is a time on the forward branch, t + is a time on the back- 
ward branch and r is a time on the imaginary branch. This 
contour is parametrized by the contour-time s, that runs from 
to t exp on the forward contour, from t exp to 2t exp on the 
backward contour and from 2t exp to 2t exp +/3 on the imaginary 
branch. Note that this mapping of imaginary- and real-time 
into the complex plane leads imaginary time (the one familiar 
from thermodynamic perturbation theory) to end up on the 
imaginary axis, therefore it's a purely imaginary number. 



With these notations we can deduce the weight of a con- 



figuration from the partition function Eq. (40) 

'-iU' "' 



W(C n ) = 



J{&et{M°(C n ))F(C n ) (44) 



where F(C n ) collects the contribution from all phases in 
the configuration: 



ncn) = n 



dz{s) 



fc=0 



ds 



(45) 



To evaluate the sum over all configurations stochastically 
we can use a Markov process using the moves of adding 
and removing a vertex. To write down their acceptance 
ratios we need, additionally to the weights, the proposal 
probabilities of the moves. The addition of a vertex is 
proposed with Tc n ^c n +i — 2WZ> which corresponds to 
the selection of an Ising spin (from {±1}), the choice of 
a site (from N sites) and of a contour-time in the range 
from [0, 2t exp + 0\. The proposal probability to remove a 
vertex is Tc n+1 ->c n — ^zp[ which corresponds to the se- 
lection of a vertex from CVi+i which has n + 1 vertices. 
We still face one issue before we can write down the ac- 
ceptance ratios for the Metropolis algorithm. Since G° in 
the real-time setting is an arbitrarily complex value and 
the expressions for the weights have explicit imaginary 
units a probabilistic interpretation is inhibited. Conse- 
quently we use their absolute values |VK(C n )| instead of 
W(C n ), but we have to compensate for this when mea- 
suring observables by keeping track of the phase of a 
configuration. We write down the acceptance ratios of 
the moves with the imaginary units still intact, keeping 



in mind that while implementing them we have to use 
the absolute values: 

/— iUNLF(C n + l ) ridct(M <T (C„ + i)) 

Pc n ^c n+1 = mm I („ + i) F (c„)ndct(A/„(c„)) > 1 

(46) 

and 

/ (n+l)F(C„)ndet(M„(C n )) 

Pc n+1 ^c n = mm I i[/JVLF ( Cri+1 )' T ndct ( MtT ( Cn+1 )), 1 

(47) 

These two moves are usually sufficient for the ergodicity 
of the algorithm. See Ref. [7] for a discussion of the cases 
where this does not apply. 



B. Measurement of observables in CTQMC 

Having generated the Markov chain of configurations 
we can start to measure observables, e.g. the single parti- 
cle Green's functions. The Green's function can be mea- 
sured by inserting the Green's function operator into the 
average Eq. (34) therefore it is necessary to evaluate the 



following expression while keeping track of the sign. 
G l3 {s,s ) = — 2^ / dz r ■■dz n x 

(TcHuiz,) . . . iJ t/ (z n )4(z( s ))c,(z( s ')))o 
E(-f)"i ? (C n )ndet(M (7 (C n ))((G ii ( S , S ')))c„ 



E 



F(C n )Y\det(M„(C n )) 



E w(c n ) 



(48) 



where we have similarly to Ref. 7\ introduced the con- 
tribution of one configuration to the observable 



„ r ( (T c H u (z 1 )...H u (z n )4(z(s))c j (z(s')))o 
{{Gij{s,s )))c = /nr tt , — s TT1 — n • 

(49) 

Now we are at the right spot to elaborate a bit on the 
sign problem. As stated before we have to replace the 
true weight W(C n ) by its absolute value |W(C n )|. We 
can repair this by rewriting the last line of |Eq. (48j with 
W(C n ) — \ W(C n )\ii(C n ). We have introduced the phase- 



factor tt(C„) = e la ^(w(c n )) 



W(C„) 



Then 



Gij(s, s ) — 



\W(Cn)\ 

EW(C n ){{G i:i (s,s')}}c n 

E w(c n ) 

c n 

E\W(Cn)\n(C n )((G i:i (s,s')))c n 

J2\w(c n )\Tr(c n ) 



(50) 
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Expanding this fraction by yj 



mc n )\ 



gives: 



E |W(C n )|7r(C„)«Gy( S , S ')» c „ 



E l^(c„)| 

On 



E |w(c„)K(c„) 
"e |w(c„)| 



(51) 



(n) ■ 

That way we see that measuring physical observables re- 
quires keeping track of the phase-afflicted observable and 
of the phase itself. The average value of the true phys- 
ical observable is then determined as their ratio. For 
the reduction of higher Green's functions to single parti- 
cle Green's functions see the algebraic identities given in 
Ref. [TB]. Equivalent methods have been published in- 
dependently |171 118] . For an overview of the algorithms 
the reader is referred to the review article jT5] ■ 



quantum Monte Carlo method provided that the stochas- 
tic evaluation of the thermal density matrix does not suf- 
fer from the negative sign problem. The starting point is 
the Trotter decomposition of the imaginary-time propa- 
gation 



M 



e -l3(H a +H v ) = nm TT e -ATH„/2 e -ArHu e -ArH 



/2 



(52) 

(At = f3/M) followed by the Hubbard-Stratonovitch de- 
coupling of the Hubbard interaction, 



-AtHu _ e 



NAtU/4 



2 N 

Sl-"SJV = ±1 



s;(n;, t -n ij4 .) 



(53) 



C. Auxiliary field quantum Monte Carlo approach 

The time evolution with a single particle Hamiltonian 
can be very efficiently computed with the auxiliary field 



with cosh(a) = e ArU ^ 2 . This allows to express the ther- 
mal density matrix - or imaginary-time propagation - in 
terms of a sum of non-interacting problems in an external 
field, at the expense of introducing a systematic error. In 
particular, for a given number of Trotter slices M, 



o N0U/4 M 



TJ e -Arff /2 e aE s (™i,T ~ n ',l) e ~ ArH /2 _|_£)(At 2 ). 



(54) 



=t/.(/9,0) 



Note that the Hubbard Stratonovitch field has acquired 
an extra imaginary-time index k. 

Within this approach the real-time expectation value 
of an observable O after a quench at time t = to a 
non- interacting Hamiltonian H reads: 



(0)(t) = 



J2 s Ti-[u s {p 7 oy tH °Oe 



itHn 



E s Tr[(7 s (/?,0)] 
_ E s Tr[^0g,O)]((Q)) s ft) 
£ s Tr [U 8 {(3, 0)] 



(55) 



with 



«°»'W = Tr[U s WM • (56) 

For a fixed Hubbard Stratonovitch configuration, Wick's 
theorem applies such that the knowledge of the single 
particle Green's function, 



[Gs(t)] x>y = (te4».(t), 



(57) 



suffices to compute the time evolution of any quantity. 
For 



H = J24lKo] x , v > 



(58) 



the Green's function matrix satisfies the equation of mo- 
tion: 



dt 



G 8 (t) = -i{H ,G s (t)}, 



such that 



G a (t) 



-itUi 



G s (t = 0)e 



itU 



(59) 



(60) 



with the Hamiltonian matrix Hq. The above equation 
reveals how to generalize standard finite temperature im- 
plementations of the auxiliary field algorithm to account 
for quenches to non-interacting Hamiltonians. In particu- 
lar for each realization of the Hubbard Stratonovitch field 
the equal time Green's function matrix, which correspond 
to the central quantity in the algorithm, can be propa- 
gated according to the above equation at the expense of 
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The double occupancy (n^n^) for different values of U 

jV = 64,/3 = 10 



0.25 



0.2 



4 



T 



T 




U = 1 

C/ = 2 

t/ = 3 

U = 4 

■ non-interacting case 
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FIG. 5. The double occupancy of the ID Monte Carlo simu- 
lations seems to decay to 0.25, the value of a free model. 



two matrix multiplications. Using Wicks theorem, arbi- 
trary correlation functions at a given time t can be com- 
puted. For an introduction to the usual imaginary-time 
formulation of the auxiliary field QMC see Ref. [8]. 



APPLICATION TO QUENCHED HUBBARD 
MODELS 



Thermalization towards a free model 



In |Fig. 5| we show the time- resolved double occupancy. 
Starting from its initial value in the Mott-insulating state 
the double occupancy shoots up to a value larger than 
that of free electrons ((n^n^) — 0.25) and peaks at a time 
of t 0.61 independently of the initially chosen U . This 
coincides nicely with the first zero of Jo (4a;). Further- 
more, the period of the following oscillations is indepen- 
dent on U . This confirms our approximate analytic result 
that the frequency of the oscillations mostly depends on 
the band width. In the long time limit (n^n^) approaches 
the non-interacting value. (Note that we limited the plot 
to a maximal time of t = 16 as this is the time scale 
where the finite size effects due to the boundary set in.) 

Since the double occupancy equilibrates to the value 
of free electrons we can conjecture that the long time 
stationary behaviour is described by an effective non- 
interacting model supplemented with Lagrange multipli- 
ers that enforce the conservation laws. We therefore pro- 



FIG. 6. Spin spin correlations in ID and 2D 



S(q = 7T,i) 



\S(q,t) - S cS (q)\,q = 7T 




(a) The spin spin correlation 
functions decay towards the 
values given by the effective 
model (dashed line). 

S(q=(n,n) T ,t) 



(b) In the log-log-plot we see 
the power-law like behaviour of 

S(ir,t). The decay is roughly 
t~ l . Here N = 64 and /3 = 10. 

\S(q,t)-S eS (q)\,q = (tt,tt) t 





(c) Also in 2D we see the decay 

towards the effective model. 
Note, that in contrast to the ID 
simulations the maximum time 
is t ss 4.5 as determined from 
finite-size scaling. 



(d) In 2D, the envelope of the 
oscillations follows a long-time 

decay law <x t~ 2 . For (7 = 8 
and /3 = 2.5, this is supported 

by the dashed line which is 
<x J Q 4 (4t) 



pose an effective density matrix of the form 

PcS = Tre-^ff^ff 
where the effective Hamiltonian is of the form 



(61) 



H, 



off 



E 

ijcr 



n k „-(n k(T {t = 0)}) (62) 



with undetermined Lagrange multipliers Afc CT and an ar- 
bitrary single particle Hamiltonian set by . With this 
ansatz one can predict uniquely the long time stationary 
value of any correlation function. In particular consider 
the equal time spin spin correlations. Owing to Wick's 
theorem they are uniquely determined by the single par- 
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ticle occupations, 



Spatially resolved charge charge correlation functions for 
different times 



aa 
kp 



(63) 



Since the single particle occupation numbers are con- 
served quantities the long time behaviour of any cor- 
relation function can be uniquely determined from the 
knowledge of (rika) at £ = 0. We test this prediction 
by computing the time-dependent spin spin correlation 
functions: 



N 



' C pa' C P 



- q ,*>)(t). (64) 



aa' kp 



Fig. 6(a)] shows the behaviour of S(q = 7T, £) in ID and 
Fig. 6(c)|in 2D. In both considered dimensions, the spin 



spin correlations decay rapidly and approach the values 
determined with the effective model. Note that in ID it 
takes longer for the oscillations to fade out than in 2D. 
This is consistent with our calculations from perturbation 
theory which predict a smaller exponent in ID than in 
2D. 



B. Decay of correlation 

To extract the decay rate of the correlation functions 
we plotted in Fig. 6(b) and Fig. 6(d) the difference to the 
effective model on a log-log scale. We see that the max- 
ima of each of the oscillations can roughly be fitted by 
straight lines, thus the decay shows power-law behaviour. 
Due to the rather low linear dimension of the 2D system 
the observable time evolution is restricted by finite-size 
effects to about t = 4.5. In ID we observe a decay with 
a power law like whereas in 2D as t~ 2 . This is con- 
sistent with the previous analytic considerations. 

To conclude we observe that the system relaxes to 
a state that is well described by a fermionic gaussian 
Hamiltonian. At least for bosons Cramer et al [2U] have 
published proofs that the time-evolution of an arbitrary 
initial state under a quadratic Bose Hamiltonian - there- 
fore some kind quench dynamics - leads to local relax- 
ation towards gaussian Hamiltonians. Physically, the au- 
thors argue that this is due to the effect that every sub- 
system acts like a bath for the other, while their coupling 
is mediated by the local interactions. 



Information propagation in correlation 
functions 



To study the information propagation in the system 
we consider two particle correlation functions. Informa- 
tion propagation has already been studied for a ID Bosc- 




FIG. 7. The spatially resolved charge charge correlation: 
\{no(t)rii(t)) — {no(t)){rii(t))\. For t = we see the char- 
acteristic exponential decay of an insulator. Between t — 1.6 
and t = 3.2 we see that a characteristic front forms that is 
propagating through the lattice. The area behind this front 
seems to be metallic as evidenced by the lack of an exponen- 
tial decay. This is a lattice of 128 sites at ft = 10 with an 
initial U — 1. 



Hubbard model in Ref. 21J and for spinless fermions in 
R ef. m. 

|Fig. 7| plots the spatially resolved charge charge cor- 
relation functions as a function of time. At t = we 
observe the characteristic exponential decay of this quan- 
tity as appropriate for insulating states. As a function 
of time a characteristic horizon forms. Beyond this hori- 
zon the charge charge correlation functions retain their 
exponential decay, whereas well within the horizon time 
independent correlation functions emerge. 

To understand the nature of the decay of the charge 



correlations well within the horizon, we plot in Fig. 8 
their time evolution for fixed distances r. As apparent, 
the equilibration time grows with the distance r as well 
as with the initial value of the Hubbard interaction U. 
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Charge charge correlation functions as a function of time for 
different values of U 





t V = i 1 


— |C„,(U)I 

— |C„„(3,t)| 

— |C„„„(5,i)| 

- 10.(1)1 

■ - lft»(3)l 

■ - |Cd,(5)| 




~ 0.1 

o 






0.01 








0.001 


1,1 




FIG. 8. For a 64 site lattice we see that the time in which the 
correlation functions Ccorr (r, t) = {no(t)n r (t)) — {no(t)){n r (t)} 
equilibrate, depends on the initial conditions, the chosen U. 
But for U = 1 it seems reasonable to think of the short- 
ranged correlation functions as equilibrated to the effective 
values C c g(r). Note that the y-axis has a logarithmic scale. 



However, the stationary value is consistent with our effec- 
tive model such that well within the horizon, the charge 
charge correlation functions are given by: 



(n(r)n(0)) cS cx ^ ]T e^ r (n ka )(l - (n k+g , a )). i (m ) 

akq 



In the above, (riko-) corresponds to the single particle 
occupation number at time t — 0. Within a mean-field 
spin density wave approximation this quantity reads: 



(nk. 



1 



<j)SDW — 



1 - 



e{k) 



(66) 



Inserting this form in Eq. |(65)| yields an exponential 
decay of the charge correlations. We note that this expo- 
nential decay of the QMC data may be very well repro- 
duced by the above equations with A ss 0.075. At our 
largest time, t — 14.4 in |Fig. 8| the charge correlations 
are converged in the region r < 16 and the maximum of 
each oscillation is consistent with an exponential decay. 
In terms of the effective model, acquaint to describing the 
long time stationary state, the SDW result of Eq. (66) 
implies that: 



Tr [pcsriko] = {nua, 



SDW 



1 

\ -|- g/3effeeff(&) ' 



(67) 



The last equation defines an effective band structure as 
well as an effective temperature [55] . With this construc- 
tion, the state after the quench may be perceived as a 
metallic state at finite temperature. 

Having discussed the velocity of the information propa- 
gation we get a rather clear-cut estimate of the time scale 



at which finite-size effects set in. In our simulations on 
lattices of 128 sites the finite-size effects set in at t « 16, 
because the velocity of the information is v « 4 and due 
to the periodic boundary conditions we can effectively 
only use half of the lattice. The torus topology of the 
lattice can be observed in Fig. 9(a)] and Fig. 9(b) where 
the horizon is symmetrically expanding from the top and 
the bottom of the figure. Calabrese and Cardy [53J HI] 
have put forward the picture that this information trans- 
port happens mainly by ballistic transport of the elec- 
trons. As we have a characteristic upper limit of the 
speed of the information propagation, this limit can be 
identified with a Lieb-Robinson bound. Lieb-Robinson 
bounds are the upper limits to the group velocities of 
excitations traveling through the considered system. As 
already mentioned, they define a light cone like structure 
that gives rise to a notion of causality, since outside of 
the cone any influence of an excitation is exponentially 
suppressed. Any non-negligible information transport is 
therefore limited by this speed. To assess if we really 
fulfill this characteristic exponential suppression of in- 
formation outside the light cone we consider some spe- 
cific values of the charge charge correlation functions as 



a function of time. In Fig. 10 we see that especially the 



longer range correlation functions show an exponential 
build-up of correlation outside of the causality cone. 

Since we see a maximum velocity of information prop- 
agation as well as the exponential suppression outside of 
the causality cone, we believe to have truly found the 
Lieb-Robinson bound in the charge sector. The spin- 
sector also shows a characteristic velocity and the expo- 
nential suppression, but our data is way more noisy for 
the spin spin correlations as is visible from Fig. 9(b) The 



fact that we could observe this finite propagation of infor- 
mation is due to our ability to do lattice simulations on 
very long chains (at least in ID), in contrast to the sim- 
ulation of an effective impurity-like model as e.g. in the 
DMFT approximation. Surprisingly little work has been 
done to include the light cone into numerical approxi- 
mation schemes, although it is a characteristic feature 
of lattice models in non-relativistic quantum mechanics. 
One exception known to the authors is Ref. [25 . 



VI. SUMMARY 

We have used two QMC methods which allow to tackle 
the general problem of quenches from correlated ther- 
mal initial states to arbitrary one-particle Hamiltonians. 
Provided that the initial density matrix can be gener- 
ated without encountering a negative sign problem, the 
real-time dynamics does not suffer from a dynamical sign 
problem. This allows to access large lattices and long 
propagation times. The algorithms used are generaliza- 
tions of the weak coupling continuous time and auxil- 
iary field QMC algorithms to the Matsubara-Keldysh 
contour. |Fig. 4| As a first step we have studied the 
dynamical transition from a Mott insulating state to a 
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FIG. 9. Comparison of the charge charge and spin spin correlation functions 




2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 

t t 
(a)We see the causality cone in the charge charge correlation (b)The spin spin correlation functions also show the horizon, but 

function. |Fig. 10| corresponds to a vertical cut along the t-axis. also more noise. This plot as well as |Fig. 9(a)| is for a 128 site 

lattice at U = f and /3 = 10. 



Exponential suppression outside of the causality cone of the 
charge charge correlation 




— n = 


8 


— n = 


24 


— n = 


32 


— n = 


40 


— n = 


56 



FIG. 10. The exponential suppression of \{no(t)n,i(t)) — 
(no(t)} {rii(t))\ outside of the causality cone. The blue lines 
are exponentials meant as a guide to the eye. The indices n 
correspond to different distances in the measurement of the 
correlation function. 



Fermi liquid corresponding to the quench from a finite 
to a vanishing Hubbard repulsion U, both in one and 
two dimensions. We find that spin spin correlation func- 
tions and double occupancy decay towards values that 
can be reproduced by an effective single particle Hamil- 
tonian where the particle densities rik a are restricted to 
the values of the thermal density matrix. We observe in 
|Fig. 6| that the decay of the magnetic order depends on 
the dimensionality of the system. In ID we observe a de- 
cay where the oscillations are enveloped by a decay like 
i -1 and in 2D as t~ 2 . Monte-Carlo methods have proven 
to be reliable tools for tackling this problem, since espe- 
cially in 2D there is a lack of non-equilibrium approxima- 
tion schemes that could provide insight into the system 



at hand. One exception known to the authors is an ex- 
tension of Cluster Perturbation theory in Ref. We 
have compared successfully our results with mean-field 
and perturbative calculations in both one and two di- 
mensions. The very good agreement points to the fact 
that the quench pumps enough energy in the isolated sys- 
tem such that the detailed correlation induced properties 
of the initial system do not effect in any significant way 
the evolution to the stationary state. This is particularly 
striking in the one-dimensional case, since the mean-field 
approximation captures by no means the physics of the 
initial Mott insulating state. Due to this large amount 
of energy released by the quench one can argue that the 
isolated system goes to a high temperature state where 
vertex corrections can be neglected. Hence any n-point 
correlation function can be described by a product of 
single particle Green's functions. Since in D dimensions 
the single particle Green's function of a non-interacting 
system exhibits a diffusive envelope, i - " 5 ", it follows that 
an n-point correlation function has a long time behav- 
ior ex t 2^ . This is confirmed by the QMC simulations 
both in one and two dimensions. In the charge charge 
correlation functions, Fig. 9(a)] we observe that the in- 
formation propagates with a velocity of v ~ 4 through 
the lattice. This behaviour is the same as predicted by 
Lieb- Robinson theorems for various systems. Thus in the 
charge charge correlation function, this system preserves 
a sense of locality. The same applies for the spin spin 
correlations, Fig. 9(b) For distances within the light 



cone, the charge charge correlations are consistent with 
a power-law decay. Beyond this length scale they follow 
an exponential law characteristic of the insulating state. 
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Appendix A: The mean-field magnetization 

1. The magnetization in ID 

The magnetization m z (A,t) in the thermodynamic 
limit reads in ID: 



m z {A,t) 



2A 

7T 



dx 



1 



1 



\/4 - x 2 a/A 2 + x 2 



cos(2ta). (Al) 



Setting t = we get the initial value of the magnetization 
in ID as 



m z (A,t = 0) = 2 K(~^] 

7T A Z 



(A2) 



where K denotes the complete elliptic integral of the first 
kind. Inserting the familiar expansion of the cosine into 
Eq. (Al) and integrating term-wise we get 



k=0 

oo 



(2fc)!fc 



2^ 2^1 



(fc!) 



H 2 



I' 

h+ i ' A 2 



A 2 



(A3) 



where («)& denotes the Pochhammer symbol and 2 ^i is 
the Gauss hypergeometric function. Using the Pfaffian 
transformation for the hypergeometric function we get 



m z {A,t) 



4 

A 2 



— i OO 



fe=0 



(fc 



n ' A 2 



(A4) 

Decomposing the Gauss hypergeometric function and us- 
ing some properties of the Pochhammer symbol we get 
again a series of hypergeometric type: 



m z (A,t) = (1 



4 ' 

A 2 , 



4 

A 2 



)" l S 2 (i,i,l,r,,-4t 2 ) (A6) 

with the definition 77 = 4(4 + A 2 ) -1 . S 2 was introduced 
by P. Humbert to denote the twice confluent version of 
Appell's F3 double hypergeometric function. See Ref. 
[27] chapter 7.2.4 for the definitions of hypergeometric 
functions in several variables. We note the expansion of 
this S2 in terms of Bessel functions: 



S 2 (^,l,r,-4t 2 ) = £ 
3=0 



(!Mik Jj(4<) . ( A7) 



2t 



3- 



So far, this particular 2 2 has resisted all attempts to 
deduce a closed form expression. Nevertheless it gives 



Comparing the asymptotic expansion to numerical 
mean-field data 




FIG. 11. A comparison of numerically gained mean-field data 
of a finite chain of 4096 sites with the asymptotic expansion 
Eq. (A8)| Here U = 2 which gives a self-consistently 



given m 
determined A 



0.34. 



the asymptotic expansion with respect to t. The right 
hand side of Eq. (A7) provides a generalized asymptotic 
series which is asymptotic for t — > 00 with respect to the 

asymptotic scale {4>j} = {t^^ 2 },j = 0, 1, Note that 

for large t, S 2 gets insensitive to changes in 77. This is due 
to the fact that the leading order behaviour of S 2 is just 
Jo(4t) without any ^-dependent prefactor. To conclude, 
we give the leading order behaviour of m z (A, t) for large 
t, 



m z (A,t) = (l + ^Y 



;j (4t). 



(A8) 



This function is plotted in |Fig. 11| for U — 2 which gives 
A 0.34. Obvious are the deviations for small t, but it is 
remarkable that the amplitude of the long time behaviour 
is very ac curately described by the inverse square root in 
Eq. (A8) 



2. The magnetization in arbitrary dimensions in 
the large A limit 

We consider the mean-field magnetization on a hyper- 
cubic isotropic lattice of dimension D. 



00 


deg D (e) 


v/ 




— 00 


OO 


2A 




— Re 


J deg D 


7T 


-00 



cos(2ie) 
Ve 2 + A 2 



(A9) 



Ve 2 + A 2 



so essentially it is just the Fourier transform of some more 
complicated function. But having rewritten it that way, 
we see that the following results also apply to the spin 
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spin correlation function as derived in perturbation the- the strip of analyticitj&is given by —1 < Re(s) < 0. Then 
ory 



Eq. (33) We start by inserting the representation . „x A f e~ 

— n J m z (t,A,D) = — / deg D {e)- 



ot 7T J Ve 2 + A 2 

— OO 

e+ioc oo 

C+ZOO 1 /" 1 /* 

A /".._„ 1 T + S, = _ h e AST{--\r(-±t\ / Jen^Ae™*** 



Vx 2 + A 2 



/d.A-^*T(-|)r ( i±^) (A10) = ^ jds^T { --n- + -) j d e 9D{ e)e^e 

V c—ioo — oo 



(All) 



in terms of a Mellin-Barnes integral (for properties of the Inserting the definition of the density of states this 
Mellin transform the reader is referred to Ref. [28 ) where lengthens to 



c+zoo 



m z (t,A,D) = ^ JdsA- s T(^ S -)T(^+ S -) J de J dk D S(e - 2 ]T C0B(Aj))e 2 "V. (A12) 



Substituting Xi — cos(ki) we get 

c+200 



c-t-ioo oo D / \ D 

m z (t,A,D) = ^ fdsA-^-nU 3 -) fdef dx D m^^-) 8(e - 2^ 

2**j 2 2 2 j J** tAW 1 -^; ti 



xAe 2Ue e s 
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(A13) 



Displacement of the integration path to the right yields 



, * n i >^r(i + fc) 



9{l-xj)\ 4itE: 
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i ^ r(| + fc) rf 2fc 

~ ttI ^ fc!4 fc cft 2fe 
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(A14) 



r 



which gives the final result 



with bandwidth w the relevant expression to analyze is 



w ah 



^'^"""^^f^^- (A15) 



m z (w,A,t) 



A 



dx 



cos(2te) 



wir J VA 2 + x 2 



(A17) 



where we have introduced the cut-off index M to tcrmi- We again perf orm an asy mptotic expansion with respect 
nate the asymptotic series. The A -> oo limit can also to A by using |Eq. (A10)| Then 
be inserted into the high-temperature expression of the 
spin spin correlation function. 



m z (w, A,t) 



1 

2W7T2 



c+ioo 



dsT(-~)T(-^-) j dx cos(2tx)x~ 



3. The magnetization in the large A limit for a 
constant density of states 

Assuming a constant density of states 



<?(e) = — 6(e 2 -«; 2 ) 
2w 



c+?oo 

f , s s ,1 + s, w 1+s A- s 
" ( -2 )r| ^»2^l(M ,f ' 



2 + 2 . j.22 



(A18) 



The poles in the positive complex half-plane are at zero 
(A16) anc ] a t even integers. Therefore displacement over the 
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.2 2 

: —t w 



first M poles at Sk — 2k yields: 

m 2 ( W ,A,0oc-^ fc!(i + 2/fc)A2fcl f 2 ^ 4 

(A19) 

The first term is the contribution at infinity, therefore 
the A — > oo limit. We get 

/ « x sm(2wt) , , 

mJw, A oo, t) = — - '-. (A20 

We proceed to evaluate the hypergeometric function by 
using an integral representation ( 27J Eq. 7.2.3.9). The 
hypergeometric function in Eq. (A19) is related to the 



spherical Bessel functions and therefore a reduction to 
a finite sum of elementary functions is possible. Setting 
x = wt we get 



77 / ^ . ,.- 1 _ r ^2 + k ) 



2r(| + fc) 



dz z cos(2irz) 



(A21) 



where we have performed a variable substitution in the 
last line and reduced the hypergeometric function in the 
integrand to a cosine. Performing integration by parts in 
the integral 2k times we get 



dz z 2k cos{2wtz) = ^^(-(2wt) 2 y k sin{2wt) 



2k-l 



7TJ, 



— ( 2k - 3 + l)i(-2^*)" j sm(2wf - t> 

W 3=0 



(A22) 



From this expression we immediately see that the only 
occurring frequency is 2w and that the slowest decay is 



for all valid values of k. 



reduced the hypergeometric function in Eq. (A19) 



As a by-product we have 
to a 



finite series of elementary functions and therefore see, 
that the oscillations are due to pure sines. 



